library(foreign)
library(RCurl)
library(xtable)


# consider changing working Directory to find this file
mydf = read.csv("data\merged_cces.csv")

sfix = function(x) {
  if(x == 2016){
    return(1)
  }
  else {
    return(0)
  }
}

mydf$treated = apply(mydf['Year'], 1,sfix)
mydf$DiD = mydf$protective*mydf$treated
mydf$DiD2 = mydf$police*mydf$treated

# import robust function from github repository
## this gives heteroskedastic "robust" standard errors equal to Stata's reg, robust
url_robust <- "https://raw.githubusercontent.com/IsidoreBeautrelet/economictheoryblog/master/robust_summary.R"
eval(parse(text = getURL(url_robust, ssl.verifypeer = FALSE)),
     envir=.GlobalEnv)


# for each variable listed above: find sample proportion by group and calculate p-values for differences
varlist = c("contacted", "VoteGOP", "defvote",  "sign" , "volunteer" , "Matched", "contribute" )
did = matrix(nrow = 7, ncol = 6, dimnames = list(varlist, c("Pop Mean (12)", "Police Mean (12)",
                                                            "Pop Mean (16)", "Police Mean (16)",
                                                            "DiD", "p-value")))

# Table 2: Political Behavior among Police Officers
idx = 1
for(var in varlist) {
  did[idx, 1] = mean(mydf[mydf$police==0 & mydf$treated==0,var])
  did[idx, 2] = mean(mydf[mydf$police==1 & mydf$treated==0,var])
  did[idx, 3] = mean(mydf[mydf$police==0 & mydf$treated==1,var])
  did[idx, 4] = mean(mydf[mydf$police==1 & mydf$treated==1,var])
  did[idx,5] = summary(lm(mydf[,var] ~ treated+police+DiD2, data=mydf), robust=TRUE)$coefficients[4,1]
  did[idx,6] = summary(lm(mydf[,var] ~ treated+police+DiD2, data=mydf), robust=TRUE)$coefficients[4,4]
  idx = idx+1
}
print(did)
# xtable(did)

# Table S2: sensitivity check
idx = 1
for(var in varlist) {
  did[idx, 1] = mean(mydf[mydf$protective==0 & mydf$treated==0,var])
  did[idx, 2] = mean(mydf[mydf$protective==1 & mydf$treated==0,var])
  did[idx, 3] = mean(mydf[mydf$protective==0 & mydf$treated==1,var])
  did[idx, 4] = mean(mydf[mydf$protective==1 & mydf$treated==1,var])
  did[idx,5] = summary(lm(mydf[,var] ~ treated+protective+DiD, data=mydf), robust=TRUE)$coefficients[4,1]
  did[idx,6] = summary(lm(mydf[,var] ~ treated+protective+DiD, data=mydf), robust=TRUE)$coefficients[4,4]
  idx = idx+1
}
print(did)
# xtable(did)



